Maxwell Mkondiwa (m.mkondiwa@cgiar.org) and Jordan Chamberlin (j.chamberlin@cgiar.org)
1 Introduction
Amidst concerns of global agricultural productivity growth slowdown, there is an emerging consensus that crop productivity growth in some African countries has either been stagnant or slow in the past decade. This slowdown is attributed to degrading soil health and volatile weather patterns resulting in low partial productivity especially of nitrogen—hereafter nitrogen use efficiency (NUE) and falling profits. We contribute to this literature by using plot level nationally representative panel data (2008-2020) for Tanzania to examine if indeed NUEs and profits have been on a downward spiral. We combine a causal random forest model for heterogeneous treatment effects and a regional market economic surplus model to explore the economic implications of the NUE trends.
The integrated framework is presented in the figure below.
This notebook provides R code for Mkondiwa and Chamberlin (2004)(https://hdl.handle.net/10883/35304) on estimating the heterogenous treatment effects in Tanzania using non-parametric geoadditive and causal machine learning (ML) models.
# Lm List by yearbaseline_quad_lmList <-lmList(yld_~N_kgha+I(N_kgha^2)+P_kgha+I(P_kgha^2)+plotha+ncrops+hhmem+femhead+headeduc+arain_tot|year,tz_lsms_panel)coef(baseline_quad_lmList )
## with interactionswith rainfall# Lm List by yearbaseline_quad_lmList_inter <-lmList(yld_~N_kgha*arain_tot+I(N_kgha^2)*arain_tot+P_kgha+I(P_kgha^2)+plotha+ncrops+hhmem+femhead+headeduc|year,tz_lsms_panel)coef(baseline_quad_lmList_inter )
visreg2d(baseline_quad_inter_soils_2010, "N_kgha", "soc_5_15cm",ylab="Soil carbon (g per kg)",xlab="Nitrogen (kg per ha)", main="2010")
Code
visreg2d(baseline_quad_inter_soils_2012, "N_kgha", "soc_5_15cm",ylab="Soil carbon (g per kg)",xlab="Nitrogen (kg per ha)", main="2012")
Code
visreg2d(baseline_quad_inter_soils_2014, "N_kgha", "soc_5_15cm",ylab="Soil carbon (g per kg)",xlab="Nitrogen (kg per ha)", main="2014")
Code
visreg2d(baseline_quad_inter_soils_2020, "N_kgha", "soc_5_15cm",ylab="Soil carbon (g per kg)",xlab="Nitrogen (kg per ha)", main="2020")
Code
# Quadratic plateau or linear plateaulibrary(easynls)tz_lsms_panel.temp <- tz_lsms_panel[, c("yld_","N_kgha")]nls.QP <-nlsfit(tz_lsms_panel.temp, model =4)nls.QP$Model
[1] "y ~ (a + b * x + c * I(x^2)) * (x <= -0.5 * b/c) + (a + I(-b^2/(4 * c))) * (x > -0.5 * b/c)"
Code
mn=nls.QP$Parametersmn
N_kgha
coefficient a -3.909774e-01
coefficient b 1.052253e-02
coefficient c -1.090000e-06
p-value t.test for a 2.841850e-01
p-value t.test for b 0.000000e+00
p-value t.test for c 0.000000e+00
r-squared 8.650000e-02
adjusted r-squared 8.640000e-02
AIC 8.165610e+04
BIC 8.168471e+04
maximum or minimum value for y 2.499715e+01
critical point in x 4.825482e+03
Code
summary(mn)
N_kgha
Min. : -0.39
1st Qu.: 0.00
Median : 0.09
Mean :14015.95
3rd Qu.: 1225.12
Max. :81684.71
5.1.1 Site-year specific Quadratic Only response functions
The site-year specific Quadratic response function can be modeled as At level 1, we have \[ Y_{i} = a_{i} + b_{i}*N + c_{i}*N^2 + \varepsilon_{i} \] At level 2, we have \[
a_{i} \sim N(a_0, \sigma_{a}^2) \\
b_{i} \sim N(b_0, \sigma_{b}^2) \\
c_{i} \sim N(c_0, \sigma_{c}^2) \\
\] This model can be estimated with the linear mixed model function lmer in R package lme4
Linear mixed model fit by REML ['lmerMod']
Formula: yld_ ~ 1 + N_kgha + N2 + (1 | year) + (0 + N_kgha | year) + (0 +
N2 | year)
Data: tz_lsms_panel.temp
REML criterion at convergence: 152692.7
Random effects:
Groups Name Std.Dev.
year (Intercept) 8.001e+01
year.1 N_kgha 2.748e+00
year.2 N2 6.472e-03
Residual 7.960e+02
Number of obs: 9426, groups: year, 5
Fixed Effects:
(Intercept) N_kgha N2
748.29010 11.24225 0.01232
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
Code
summary(lmer.Q)
Linear mixed model fit by REML ['lmerMod']
Formula: yld_ ~ 1 + N_kgha + N2 + (1 | year) + (0 + N_kgha | year) + (0 +
N2 | year)
Data: tz_lsms_panel.temp
REML criterion at convergence: 152692.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.3333 -0.6213 -0.3108 0.2772 6.6397
Random effects:
Groups Name Variance Std.Dev.
year (Intercept) 6.402e+03 8.001e+01
year.1 N_kgha 7.554e+00 2.748e+00
year.2 N2 4.188e-05 6.472e-03
Residual 6.337e+05 7.960e+02
Number of obs: 9426, groups: year, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 748.29010 36.88270 20.288
N_kgha 11.24225 1.75175 6.418
N2 0.01232 0.01525 0.808
Correlation of Fixed Effects:
(Intr) N_kgha
N_kgha -0.051
N2 0.046 -0.657
fit warnings:
Some predictor variables are on very different scales: consider rescaling
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Or, althernatively, can be estimated with the non-linear mixed model function nlme
Code
library(nlme)nlme.Q <-nlme(yld_ ~ (a + b * N_kgha + c * (N_kgha^2)),data = tz_lsms_panel,fixed = a + b + c ~1,random = a + b + c ~1,groups =~ year,start =c(800, 10, -0.001))nlme.Q
Nonlinear mixed-effects model fit by maximum likelihood
Model: yld_ ~ (a + b * N_kgha + c * (N_kgha^2))
Data: tz_lsms_panel
Log-likelihood: -76345.33
Fixed: a + b + c ~ 1
a b c
747.82829275 11.16211183 0.01323599
Random effects:
Formula: list(a ~ 1, b ~ 1, c ~ 1)
Level: year
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
a 57.45976820 a b
b 1.64304646 -0.022
c 0.02343941 0.946 -0.341
Residual 795.98608057
Number of Observations: 9426
Number of Groups: 5
5.1.2 Site-year specific Quadratic-plus-plateau response functions
The site-year specific response function can be modeled using a hierarchical quadratic-plus-plateau model. At level 1, we have \[ Y_{i} = \min(a_{i} + b_{i}*N + c_{i}*N^2, Y_{max}) + \varepsilon_{i} \] At level 2, we have \[
a_{i} \sim N(a_0, \sigma_{a}^2) \\
b_{i} \sim N(b_0, \sigma_{b}^2) \\
c_{i} \sim N(c_0, \sigma_{c}^2) \\
Y_{max} = a_{i} - b_{i}^2/(4*c_i)
\]
It seems the R function nlme would work to estimate this model
Code
# Define quadratic-plus-plateau functionmq <-lm(yld_ ~ N_kgha +I(N_kgha^2), data=tz_lsms_panel)a0 <-coef(mq)[[1]]b0 <-coef(mq)[[2]]c0 <-coef(mq)[[3]]clx0 <--0.5*b0/c0# Test nls# fx.QP <- function(N, a, b, c) {# y <- (a + b * N + c * I(N^2)) * (N <= -0.5 * b/c) + (a + I(-b^2/(4 * c))) * (N > -0.5 * b/c)# return(y)# }# # nls.QP <- nls(Y ~ fx.QP(N, a, b, c),# start = list(a = a0, b = b0, c = c0), data = dat.Puntel.CC.mean,# control = nls.control(maxiter = 6000))# quadplat <- function(x, a, b, clx) {# ifelse(x < clx, a + b * x + (-0.5*b/clx) * x * x, # a + b * clx + (-0.5*b/clx) * clx * clx)# }# # nls.QP <- nls(Y ~ quadplat(N, a, b, clx), # start = list(a = a0, b = b0, clx = clx0), data = dat.Puntel.CC.mean, # control = nls.control(maxiter = 6000))# nlme.QP <- nlme(Y ~ fx.QP(N, a, b, c),# data = dat.Puntel.CC.mean,# fixed = a + b + c ~ 1,# random = a + b + c ~ 1,# groups = ~ year,# start = c(a0, b0, c0))# # nlme.QP# tz_lsms_panel.nlme <- nlme(yld_ ~ (a + b * N_kgha + c * I(N_kgha ^2)) * (N_kgha <= (-0.5 * b/c)) + (a- I(b^2/(4 * c))) * (N_kgha > (-0.5 * b/c)),# data = tz_lsms_panel,# fixed = a + b + c ~ 1,# random = a + b + c ~ 1,# groups = ~ year,# start = c(a0, b0, c0))# tz_lsms_panel.nlme
Best linear fit using forest predictions (on held-out data)
as well as the mean forest prediction as regressors, along
with one-sided heteroskedasticity-robust (HC3) SEs:
Estimate Std. Error t value Pr(>t)
mean.forest.prediction 0.99724 0.10227 9.7509 < 2.2e-16 ***
differential.forest.prediction 1.25318 0.26830 4.6708 1.521e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Best linear fit using forest predictions (on held-out data)
as well as the mean forest prediction as regressors, along
with one-sided heteroskedasticity-robust (HC3) SEs:
Estimate Std. Error t value Pr(>t)
mean.forest.prediction 0.99353 0.06987 14.2196 < 2.2e-16 ***
differential.forest.prediction 1.23712 0.26364 4.6925 1.369e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
library(ggridges)library(dplyr)library(ggplot2)tau.multi_fert.forest_continuous <-predict(multi_fert.forest_continuous, target.sample ="all", estimate.variance =TRUE)tau.multi_fert.forest_continuous <-as.data.frame(tau.multi_fert.forest_continuous)tau.multi_fert.forest_X <-data.frame(tz_lsms_panel_estim_fert, tau.multi_fert.forest_continuous)NUE_Density=ggplot(tau.multi_fert.forest_X, aes(x = predictions, y ="", fill =factor(stat(quantile)))) +stat_density_ridges(geom ="density_ridges_gradient", calc_ecdf =TRUE,quantiles =4, quantile_lines =TRUE ) +scale_y_discrete(expand =c(0.01, 0)) +scale_fill_viridis_d(name ="Quartiles") +expand_limits(y =1) +theme_bw(base_size =16) +labs(x ="N use effect", y ="Density")NUE_Density
7 Exploring heterogeneity in Conditional N Use Efficiencies from CRF
7.1 Causal RF
Code
library(ggplot2)NperHa_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$N_kgha >0), aes(N_kgha, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Applied N per ha", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))NperHa_CATE_N_plot
Code
Plotsize_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$plotha >0), aes(plotha, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Plot size (ha)", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))Plotsize_CATE_N_plot
Code
P_CATE_N_plot <-ggplot(tau.multi_fert.forest_X, aes(P_kgha, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="P", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))P_CATE_N_plot
Code
Hhmem_CATE_N_plot <-ggplot(tau.multi_fert.forest_X, aes(hhmem, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="HHMEM", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))Hhmem_CATE_N_plot
Code
# By yeartau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==1]=2008tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==2]=2010tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==3]=2012tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==4]=2014tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==5]=2020tau.multi_fert.forest_X$year_det=as.numeric(tau.multi_fert.forest_X$year_det)NperHa_CATE_N_plot_yr <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$N_kgha >0), aes(N_kgha, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +lims(x =c(0, 100)) +labs(x ="Applied N per ha", y ="N use efficiency") +theme_bw(base_size =16) +facet_wrap(~year_det)NperHa_CATE_N_plot_yr
Code
# By soil organic carbonlibrary(ggplot2)soc_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$soc_5_15cm >0), aes(soc_5_15cm, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Soil organic carbon (g per kg)", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))soc_CATE_N_plot
Code
ggsave("figures/soc_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# By soil sandlibrary(ggplot2)sand_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$sand_0_5cm >0), aes(sand_0_5cm, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Sand (%)", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))sand_CATE_N_plot
Code
ggsave("figures/sand_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# Electrical conductivitylibrary(ggplot2)soil_elec_cond_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$ECN_5_15cm >0), aes(ECN_5_15cm, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Soil electrical conductivity", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))soil_elec_cond_CATE_N_plot
Code
ggsave("figures/soil_elec_cond_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# pHlibrary(ggplot2)soil_pH_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$pH_0_5cm >0), aes(pH_0_5cm, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Soil pH", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))soil_pH_CATE_N_plot
Code
ggsave("figures/soil_pH_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# soil nitrogenlibrary(ggplot2)soil_nitrogen_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$nitrogen_0_5cm >0), aes(nitrogen_0_5cm, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Soil nitrogen (g per kg)", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))soil_nitrogen_CATE_N_plot
Code
ggsave("figures/soil_nitrogen_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# By population densitylibrary(ggplot2)pop_density_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$population_density >0), aes(population_density, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Population density", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))pop_density_CATE_N_plot
Code
ggsave("figures/pop_density_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# By elevationlibrary(ggplot2)elev_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$wc2.1_30s_elev >0), aes(wc2.1_30s_elev, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Elevation (masl)", y ="N use efficiency")previous_theme <-theme_set(theme_bw(base_size =16))elev_CATE_N_plot
Code
ggsave("figures/elev_CATE_N_plot.png", dpi=300, width=4.16,height=4.16)# By rainfalllibrary(ggplot2)rainfall_CATE_N_plot <-ggplot(subset(tau.multi_fert.forest_X, tau.multi_fert.forest_X$arain_tot >0), aes(arain_tot, predictions)) +geom_smooth(method ="loess", formula = y ~ x, col ="darkblue") +labs(x ="Rainfall (mm)", y ="N use efficiency")+facet_wrap(~year_det)previous_theme <-theme_set(theme_bw(base_size =16))rainfall_CATE_N_plot
8 Are N use efficiencies falling overtime?
8.1 Causal RF
Code
tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==1]=2008tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==2]=2010tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==3]=2012tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==4]=2014tau.multi_fert.forest_X$year_det[tau.multi_fert.forest_X$year==5]=2020tau.multi_fert.forest_X$year_det=as.numeric(tau.multi_fert.forest_X$year_det)Year_CATE_N_plot=ggplot(tau.multi_fert.forest_X,aes(year_det,predictions))+geom_smooth(method="loess",formula=y~x,col="darkblue")+labs(x="Year",y="N use efficiency")+scale_x_continuous(breaks =c(2008, 2010, 2012,2014,2020))previous_theme <-theme_set(theme_bw())Year_CATE_N_plot
library(ggridges)library(dplyr)library(ggplot2)tau.multi_fert.forest_dummy <-predict(multi_fert.forest_binary, target.sample ="all", estimate.variance =TRUE)tau.multi_fert.forest_dummy <-as.data.frame(tau.multi_fert.forest_dummy)tau.multi_fert.forest_X_dummy <-data.frame(tz_lsms_panel_estim_fert, tau.multi_fert.forest_dummy)ggplot(tau.multi_fert.forest_X_dummy, aes(x = predictions, y ="", fill =factor(stat(quantile)))) +stat_density_ridges(geom ="density_ridges_gradient", calc_ecdf =TRUE,quantiles =4, quantile_lines =TRUE ) +scale_y_discrete(expand =c(0.01, 0)) +scale_fill_viridis_d(name ="Quartiles") +expand_limits(y =1) +theme_bw(base_size =16) +labs(x ="Inorgatic fert use effect", y ="Density")
Code
#library(sp)library(sf)tau.multi_fert.forest_X_dummy_sp <-SpatialPointsDataFrame(cbind(tau.multi_fert.forest_X_dummy$lon_modified, tau.multi_fert.forest_X_dummy$lat_modified), data = tau.multi_fert.forest_X_dummy, proj4string =CRS("+proj=longlat +datum=WGS84"))library(tmap)tmap_mode("plot")Nuse_effect_map <-tm_shape(tau.multi_fert.forest_X_dummy_sp) +tm_dots(col ="predictions", title ="Effect of N use on yield (kg/ha)", style ="quantile") +tm_layout(legend.outside =TRUE)Nuse_effect_map
Code
tmap_save(Nuse_effect_map, "figures/Nuse_effect_map.png", width =600, height =600, asp =0)# N Use efficiency maplibrary(sp)library(sf)tau.multi_fert.forest_X_sp <-SpatialPointsDataFrame(cbind(tau.multi_fert.forest_X$lon_modified, tau.multi_fert.forest_X$lat_modified), data = tau.multi_fert.forest_X, proj4string =CRS("+proj=longlat +datum=WGS84"))library(tmap)tmap_mode("plot")Nuse_efficiency_map <-tm_shape(tau.multi_fert.forest_X_sp) +tm_dots(col ="predictions", title ="NUE (kg Maize/Kg N)", style ="quantile") +tm_layout(legend.outside =TRUE)Nuse_efficiency_map
Code
tmap_save(Nuse_efficiency_map, "figures/Nuse_efficiency_map.png", width =600, height =600, asp =0)# N Use maplibrary(sp)library(sf)tau.multi_fert.forest_X_sp_small <-subset(tau.multi_fert.forest_X_sp, tau.multi_fert.forest_X_sp$N_kgha >0)library(tmap)tmap_mode("plot")Nuse_map <-tm_shape(tau.multi_fert.forest_X_sp_small) +tm_dots(col ="N_kgha", title ="N applied (Kg/ha)", style ="quantile") +tm_layout(legend.outside =TRUE)Nuse_map
Code
tmap_save(Nuse_map, "figures/Nuse_map.png", width =600, height =600, asp =0)
tau.multi_fert.forest_X_sp$region_name <-over(tau.multi_fert.forest_X_sp, TZ_sp[,"NAME_1"])# Summary by region and by yearlibrary(modelsummary)tau.multi_fert.forest_X_sp_dt <-tau.multi_fert.forest_X_sp@datatau.multi_fert.forest_X_sp_dt$region_name_final=tau.multi_fert.forest_X_sp_dt$region_name$NAME_1NUE_cont <-datasummary(Heading("region_name")*region_name_final~Heading("N_obs") * N +Heading("percent") *Percent()+ predictions*(Mean+SD),data = tau.multi_fert.forest_X_sp_dt, output="data.frame")NUE_cont$N_obs=as.numeric(NUE_cont$N_obs)NUE_cont$Mean=as.numeric(NUE_cont$Mean)NUE_cont$SD=as.numeric(NUE_cont$SD)library(reactable)library(htmltools)library(fontawesome)htmltools::browsable(tagList( tags$button(tagList(fontawesome::fa("download"), "Download as CSV"),onclick ="Reactable.downloadDataCSV('NUE_cont', 'NUE_cont.csv')" ),reactable( NUE_cont,searchable =TRUE,defaultPageSize =38,elementId ="NUE_cont" ) ))
Code
# Remove regions with fewer observations than 50NUE_cont=subset(NUE_cont,NUE_cont$N_obs>50)NUE_cont_sf=merge(TZ_sf,NUE_cont, by.x="NAME_1",by.y="region_name")# Map the library(tmap)tmap_mode("plot")tm_shape(NUE_cont_sf) +tm_polygons(col ="Mean", title ="Average NUE", style ="quantile") +tm_layout(legend.outside =TRUE)
pred <-SpatialPoints(elevationglobal_geodata_aoi)@coordspred <-as.data.frame(pred)names(pred)[1:2] <-c("lon", "lat")pred <- pred %>%drop_na()library(sp)pred_sp <-SpatialPointsDataFrame(cbind(pred$lon, pred$lat), data = pred, proj4string =CRS("+proj=longlat +datum=WGS84"))yield_gain_hat <-predict(b, newdata = pred)yield_gain_hat <-as.data.frame(yield_gain_hat)pred_yield_gain_hat <-cbind(pred, yield_gain_hat )pred_yield_gain_hat$sigma <-NULLlibrary(terra)yield_gain_ras <-rast(pred_yield_gain_hat, type ="xyz")plot(yield_gain_ras )
Code
library(raster)yield_gain_ras2 <-raster(yield_gain_ras)library(sf)TZ_sf_dis_sp <-as_Spatial(TZ_sf_dis)yield_gain_ras2 <-mask(yield_gain_ras2, TZ_sf_dis_sp)plot(yield_gain_ras2, main ="Yield gains to N")
